£5Llf  2*2 


PL-TR-92-2158 


AD-A258  597 


INVESTIGATIONS  OF  JOINT  SEISMIC  AND  ELECTROMAGETIC 
METHODS  FOR  NUCLEAR  TEST  MONITORING 


Charles  B.  Archambeau 
John  B.  Davies 
Jeffrey  Orrey 


University  of  Colorado 
Department  of  Physics/TAGG 
Campus  Box  583 
Boulder,  CO  80309 


February  15,  1992 


OTIC 

ELECTE 
OCT  0  81992 

A 


Scientific  Report  No.  1 


Approved  for  public  release;  distribution  unlimited 


PHILLIPS  LABORATORY 

AIR  FORCE  SYSTEMS  COMMAND 

HANSCOM  AIR  FORCE  BASE,  MASSACHUSETTS  01731-5000 


The  views  and  conclusions  contained  in  this  document  are  those  of  the  authors 
and  should  not  be  interpreted  as  representing  the  official  policies,  either 
expressed  or  implied,  of  the  Air  Force  or  the  U.S.  Government. 

This  technical  report  has  been  reviewed  and  is  approved  for  publication. 


Solid  Earth  Geophysics  Branch 
Earth  Sciences  Division 


Solid  Earth  Geophysics  Branch 
Earth  Sciencs  Division 


This  document  has  been  reviewed  by  the  ESD  Public  Affairs  Office  (PA)  and  is 
releasable  to  the  National  Technical  Information  Service  (NTIS). 

Qualified  requestors  may  obtain  additional  copies  from  the  Defense  Technical 
Information  Center.  All  others  should  apply  to  the  National  Technical 
Information  Service. 

If  your  address  has  changed,  or  if  you  wish  to  be  removed  from  the  mailing  list, 
or  if  the  addressee  is  no  longer  employed  by  your  organization,  please  notify 
PL/IMA,  Hanscom  AFB  MA  01731-5000.  This  will  assist  us  in  maintaining  a 
current  mailing  list. 

Do  not  return  copies  of  this  report  unless  contractual  obligations  or  notices  on 
a  specific  document  requires  that  it  be  returned. 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  NO.  0704-0188 


P  ,w.c  -toorjn-;  ry-'ttr.  >fr  >*>,  'O'Vclicn  of  m'jrr-f  ■  :  t.  "mir  '  -  u- 


the  io»  reviewing  •''W.ctiOF't.  serening  e«  *:-«g  c* :a  sour<*v 


aathef wig  *nd  n«tni*inifsg  th*  d»ta  needed.  »nd  cc^o*ct'«q  rf'HJ  /ev»e»v.nq  the  collection  ot  tnlo^maTiDn  Send  comment  regarding  th»*  burden  estimate  or  »n«  other  jsoect  of  this 
collection  of  information,  including  suggestions  for  re:-cmg  mu  Surom  to  Wasmnoion  Headauaners  Services.  Directorate  tor  information  Operations  and  Reoons.  12  IS  Jfffeoon 
Davis  Highway  Suite  U04  Arlington.  vA  22202-4302  and  to  me  Office  o»  Management  and  Budget.  RioerworR  Reduction  Project  (07Q4-0>881.  Washington  OC  2CS03 


1.  AGENCY  USE  ONLY  (Leave  blank) 


4.  TITLE  ANO  SUBTITLE 


2.  REPORT  OATE 

February  15,  1992 


3.  REPORT  TYPE  ANO  DATES  COVERED 

Scientific  No.  1 


Investigations  of  Joint  Seismic  and  Electromagnetic 
Methods  for  Nuclear  Test  Monitoring 


6.  AUTHOR(S) 

Charles  B.  Archambeau 
John  B.  Davies 
Jeffrey  Orrey 


7.  PERFORMING  ORGANIZATION  NAME(S)  ANO  AODRESS(ES) 

University  of  Colorado 
Dept,  of  Physics/TAGG 
Campus  Box  583 
Boulder,  CO  80309 


5.  FUNDING  NUMBERS 

Contract  No. 
F19628-9G-K-005 1 
PE  62101R 
PR  7600 
TA  09 
WUBJ 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


TAGG  001 


9.  SPONSORING  /MONITORING  AGENCY  NAME(S)  ANO  AOORESS(ES) 

Phillips  Laboratory 
Hanscom  AFB 

Massachusetts  01731-5000 

Contract  Manager:  James  Lewkowicz/GPEH 


10.  SPONSORING/ MONITORING 
AGENCY  REPORT  NUMBER 


PL-TR-92-2158 


12a.  DISTRIBUTION /AVAILABILITY  STATEMENT 


12b  DISTRIBUTION  CODE 


approved  for  public  release,  distribution  unlimited 


13.  abstract  (Maximum 200 wordi)  This  study  is  designed  to  develop  both  linear  and  non-linear 
modeling  methods  for  prediction  of  excitation  of  atmospheric  and  seismic  disturbances 
from  near  surface  explosions  and  earthquake  sources  in  complex  media  which  include 
strong  lateral  variability,  randomness  and  non-linear  response  effects.  Comparisons 
with  observations  are  systematically  pursued  to  evaluate  these  models  and  to  develop 
source  discrimination  methods.  In  this  report  we  describe  examples  of  non-linear  at¬ 
mospheric  excitation  by  near  surface  explosions  which  are  carried  to  ionospheric 
heights  in  order  to  predict  fluctuations  in  ionospheric  electron  densities  and  ioni¬ 
zation  layer  positions.  These  predictions  are  compared  to  active  EM  monitoring  by 
ground  stations  of  ionospheric  layer  perturbations  due  to  the  large  amplitude  atmos¬ 
pheric  waves  from  surface  explosions.  We  also  consider  atmospheric  turbulent  couplin; 
at  the  earth's  free  surface  and  investigate  the  high  frequency  seismic  noise  charac¬ 
teristics  resulting  from  this  source  using  combined  seismic  and  atmospheric  modeling. 
Results  are  in  first  order  agreement  with  spectral  observations  of  high  frequency 
noise,  in  that  both  the  theoretical  and  observed  particle  velocity  spectra  vary  as 
1/f  in  the  high  frequency  range  from  I  to  50  Hz.  Examples  of  the  effects  of  random¬ 
ness,  large  scale  lateral  variations  and  near  source  surface  topography  on  the  seis¬ 
mic  wave  field  from  near  surface  explosions  are  systematically  investigated.  We  find 
good  first  order  correlations  of  predictions  with  the  complex  regional  and  near 
field  observations  using  reasonable  structure/topographic  models. 


14  SUBJECT  TERMS 

Seismic  Waves 
Electromagnetic  waves 


Numerical  modeling 
Seismic  Discrimination 


15.  NUMBER  OF  PAGES 

72 


16.  PRICE  COOE 


17.  SECURITY  CLASSIFICATION  118.  SECURITY  CLASSIFICATION  119.  SECURITY  CLASSIFICATION  I  20.  LIMITATION  OF  ABSTRACT 
OF  REPORT  I  OF  THIS  PAGE  OF  ABSTRACT 


Unclassified 


NSN  7540-01-280-5500 


Unclassified 


Unclassified 


Standard  Form  298  (Rev  2-89) 

*re%crif»d  ANV  Sw 

m- 102 


Investigations  of  Joint  Seismic  and  Electromagnetic 
Methods  for  Nuclear  Test  Monitoring 


Table  of  Contents 


I.  Introduction:  Objectives  and  Methods .  1 

II.  High  Frequency  Seismic  Signal  Detection  and  Spectral  Discrimination: 


Observations  .  7 

in.  Modeling  High  Frequency  Seismic  Signal  Propagation .  18 


IV.  Modeling  Atmospheric  Wave  Fields  and  Ionospheric  Electron  Density 


Variations  Due  to  Near  Surface  Seismic  Sources  .  40 

V.  Modeling  High  Frequency  Seismic  Noise:  Atmospheric  Sources  .  52 

VI.  Summary  and  Conclusions  .  54 

References  .  59 


Acceeion  r or 

_  1 

N7!S 

OTIC 

cfmi 

TAG 

y 

"i 

U:,an- 

JUitifj 

OU:  ;oc'd 
nation 

U 

By 

Dtst.  ibution/  j 

Availability  Codes  j 

Dbt 

Avail 

$p<- 

i  ■'<  •  or 
r.i.'il 

/)-/ 

j 

iii . 


I.  Introduction:  Objectives  and  Methods 


Objectives 

This  study  is  designed  to  develop  both  linear  and  nonlinear  wave  propagation 
methods  that  can  model  the  excitation  and  propagation  of  atmospheric  and  seismic 
waves  from  explosion  and  earthquake  sources  in  realistic,  complex  media  models 
which  include  strong  lateral  variability,  randomness  and  nonlinear  response  effects. 

In  modeling  the  excitation  of  the  atmosphere  and  ionosphere  we  include  the  usual 
non-linear  transport  effects  as  will  as  ionization  in  order  to  infer  secondary  effects  pro¬ 
duced  by  large  amplitude  neutral  wave  propagation  upward  into  the  ionosphere  from  a 
surface  or  near  surface  explosion.  The  overall  objectives  of  the  atmospheric- 
ionospheric  modeling  are  (a.)  to  predict  fluctuations  in  the  electron  densities  and  ioni¬ 
zation  layer  positions  in  the  ionosphere  which  can  be  correlated  with  active  EM  moni¬ 
toring  by  ground  stations  and  (b.)  to  predict  secondary  EM  field  emissions  from  ion 
and  electron  movements  induced  by  the  large  amplitude  atmospheric  waves  from 
below.  Here  the  idea  is  to  evaluate  and  design  active  and  passive  EM  sensing 
methods  coupled  with  seismic  methods  to  define  a  monitoring  environment  which  will 
allow  large  industrial  explosions  to  be  easily  identified  based  on  the  strengths  and 
character  of  the  seismic  and  atmospheric-ionospheric  disturbances  producer". 

The  objectives  of  the  seismic  wave  propagation  modeling  are  to  take  account  of 
near  source  non-linear  phenomena  as  well  as  topographic  effects,  medium  randomness 
and  strong  lateral  variability  in  the  earth  structure,  particularly  in  the  crust  and  upper 
mantle.  We  hope  to  obtain  close  fits  to  the  complex  seismic  wave  fields  observed  at 
regional  and  teleseismic  distances  and  in  so  doing,  to  generate  a  basis  for  refined 
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detection  and  discrimination  of  small  seismic  events. 


Basic  Concepts  and  Approach:  Electromagnetic  Methods 

If  small  nuclear  tests  are  detonated  in  a  decoupling  cavity,  then  their  signals  in 
the  low  frequency  range  below  5 Hz  are  reduced  by  nearly  two  orders  of  magnitude. 
In  this  case  the  decoupled  nuclear  explosion  produces  signals  of  the  same  size  as  com¬ 
mon  industrial  explosions,  of  which  there  are  many  thousands  per  year  in  industrial 
areas.  Consequently,  it  will  be  necessary  to  be  able  to  seismically  distinguish  between 
these  numerous  industrial  explosions  and  possible  decoupled  nuclear  tests  if  a  treaty 
banning  such  tests  were  to  rely  principally  on  seismic  methods  for  verification.  At  the 
present  time  there  is  no  well  documented  method  for  such  discrimination,  although  it 
is  likely  that  seismic  (spectral)  methods  employing  new  high  frequency  (.5  to  50  Hz) 
seismic  detectors,  operating  with  very  low  internal  noise  and  deployed  at  depth  or  as 
"tight  arrays”  to  reduce  high  frequency  earth  noise,  will  make  it  possible  to  distinguish 
between  these  types  of  explosions,  as  well  as  between  small  earthquakes  and  tamped 
or  decoupled  nuclear  explosions.  In  this  study  we  seek  to  build  a  firm  understanding 
of  the  regional  seismic  wave  fields  produced  by  different  types  of  seismic  sources  in 
order  to  properly  define  discrimination  methods  and  procedures  and  be  able  to  test  and 
predict  their  variability  and  sensitivity  in  different  regional  structures. 

In  addition  to  seismic  methods  for  event  identification,  there  are  other  possibilities 
that  are  beginning  to  receive  serious  consideration.  Clearly,  sensing  of  acoustic  signals 
from  large  industrial  explosions  is  a  possible  means  of  identifying  these  events,  since 
neither  nuclear  tests  nor  earthquakes  will  produce  such  a  large  signal  in  the  atmo¬ 
sphere  However,  because  of  signal  attenuation  and  high  acoustic  noise  levels  at  the 
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earth’s  surface,  even  relatively  large  acoustic  signals  from  most  industrial  explosions 
of  interest  may  not  be  observable  beyond  a  few  hundred  kilometers.  Nevertheless, 
acoustic  sensors  located  quite  near  active  mining  areas  would  be  useful  in  identifying 
the  large  near  surface  explosions  that  are  of  greatest  importance.  Since  large  scale 
mining  areas  are  relatively  rare,  one  would  need  to  monitor  only  a  few  areas  using  this 
method  to  achieve  nearly  total  monitoring  coverage  of  the  largest  industrial  explosions. 
Consequently  locally  distributed  acoustic  sensor  arrays  around  major  mining  areas  can 
provide  critical  data  for  monitoring,  particularly  when  coupled  with  similar  seismic 
monitoring  arrays. 

Since  even  small  industrial  explosions  at  normal  mining  depths  will  produce  a 
much  larger  signal  in  the  atmosphere  than  a  decoupled  or  tamped  nuclear  test,  it  is 
possible  that  the  existance  of  this  large  signal  difference  could  be  used  in  other  ways, 
besides  direct  acoustic  field  measurents,  to  help  identify  these  sources.  That  is,  while 
even  the  larger  acoustic  signals  will  be  small  relative  to  noise  after  propagation  over  a 
few  hundred  kilometers  in  the  atmosphere  near  the  earth’s  surface,  this  is  not  the  case 
for  the  acoustic  wave  field  propagated  directly  upward  into  the  upper  atmosphere  and 
ionosphere.  In  this  case  the  fact  that  the  air  density  decreases  rapidly  with  altitude 
causes  the  amplitude  of  an  acoustic  wave  to  increase  rapidly  with  increasing  height 
and  to  strongly  perturb  the  ionosphere  over  an  area  of  the  order  of  100  kilometers  in 
radius  around  the  source  epicenter.  Consequently,  remote  EM  sensing  of  the  iono¬ 
sphere  using  radio  frequency  transmitter  -  receiver  systems  can  detect  ionospheric  per¬ 
turbations  due  to  the  acoustic  waves  from  industrial  explosions.  Observationally  the 
acoustic  waves  produce  ionospheric  boundary  motions  resulting  in  doppler  shifts  in  the 
reflected  EM  signals  recorded  on  the  ground.  By  placing  radio  frequency  receivers 
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and  transmitters  in  a  distributed  network,  comparable  to  an  in  -  country  seismic  moni¬ 
toring  network  of  about  30  stations,  it  appears  possible  to  provide  complete  monitoring 
capability  on  a  continental  scale.  In  principal  it  should  be  possible  to  identify  large 
industrial  explosions  with  high  probability  due  to  the  EM  signal  shifts  observed.  Cou¬ 
pled  with  seismic  monitoring  then,  the  occurrence  of  an  event  that  had  an  explosive 
seismic  signature,  but  that  produced  no  ionospheric  EM  effect,  would  indicate  a  prob¬ 
able  decoupled  or  tamped  nuclear  test. 

Besides  the  possibility  of  active  monitoring  of  the  strong  perturbations  in  the 
reflecting  layers  of  the  'onosphere  by  radio  frequency  sounding,  it  is  possible  to  detect 
secondary  EM  emissions  from  the  ionosphere.  Thus,  passive  electromagnetic  monitor¬ 
ing  of  the  electromagnetic  environment,  particularly  at  low  radio  frequencies,  offers  yet 
another  opportunity  to  address  source  identification  and  discrimination  issues. 

As  a  consequence  of  these  possibilities,  we  focus  on  modeling  the  atmospheric 
and  ionospheric  disturbances  produced  by  near  surface  explosions  of  various  types  in 
order  to  provide  an  understanding  and  quantitative  prediction  of  the  magnitude  and 
character  of  these  effects.  Based  at  least  partly  on  such  results,  we  can  then  hope  to 
define  and  test  particular  methods  of  acoustic  and  EM  monitoring  that  can  be  effective. 

Basic  Concepts  and  Approach:  Seismic  Methods 

The  seismic  discrimination  approach  envisioned  is  to  make  use  of  the  enriched 
high  frequency  content  of  seismic  waves  radiated  by  decoupled  and  tamped  nuclear 
explosions  as  compared  to  earthquakes  and  large  industrial  explosions.  Here  the  idea 
is  that  earthquakes,  producing  the  same  low  frequency  compressional  signal  level  as 
(say)  a  decoupled  nuclear  explosion,  will  produce  lower  levels  of  high  frequency 
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compressional  wave  energy  and  so  the  ratio  of  high  frequency  compressional  (P)  v-ave 
energy  to  low  frequency  P  wave  energy  will  be  different  for  the  two  source  types. 
Physically  this  occurs  because  the  source  dimension  of  the  comparable  earthquake  is 
much  larger  than  that  for  the  explosion  and,  further,  because  the  energy  release  process 
is  fundamentally  different.  (This  is  discussed  in  detail  by  Evemden,  Archambeau  and 
Cranswick,  Rev.  of  Geophysics,  1986). 

On  the  other  hand,  while  a  decoupled  nuclear  explosion  and  a  large  industrial 
explosion  release  energy  from  roughly  comparable  volumes,  the  nonlinear  effects 
within  the  source  volume  are  expected  to  be  sufficiently  different  to  allow  discrimina- 
non  using  compressional  wave  spectral  differences  of  the  same  sort  that  are  used  to 
discriminate  nuclear  tests  and  earthquakes. 

In  this  regard,  industrial  explosions  are  distributed  spatially,  with  partial  charges 
typically  in  bore  holes  with  regular  spacing  (tens  of  meters  intervals)  and  often 
detonated  with  time  delays  between  the  firing  of  individual  charges.  This  is  done  to 
avoid  massive  localized  surface  rupture  and  "blow-out"  and  also  to  enhance  fracturing 
of  the  material  over  the  largest  possible  volume  using  the  smallest  total  amount  of 
charge.  As  one  important  consequence,  the  time-space  separation  of  the  charges  will 
result  in  nonlinear  constructive  interference1^  of  the  shock  waves  from  the  individual 
charges.  This  interference  results  in  maximum  localized  over  -  pressures  and  very 
efficient  fracturing  of  the  material.  When  viewed  as  a  seismic  source,  that  is  as  a  radi¬ 
ator  of  far-field  elastic  waves,  such  a  source  expends  a  larger  proportion  of  it’s  total 
energy  in  fracturing  the  surrounding  material  than  would  the  same  total  charge 

t  interacting  shock  waves  will  form  "Mach  Stems"  corresponding  to  a  moving  plane  of  shock  intersection,  where  the 
pressure  jump  is  significantly  larger  than  either  of  the  pressure  jumps  along  the  individual  shock  fronts. 
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detonated  in  a  single  small  volume  within  the  medium.  Compared  to  a  tamped  nuclear 
explosion,  or  to  a  decoupled  nuclear  test,  the  amount  of  energy  radiated  as  seismic 
waves  should  be  significantly  lower  when  normalized  for  yield.  More  specifically,  the 
microfracturing  produced  should  greatly  reduce  the  high  frequency  energy  available  to 
be  radiated  from  a  spatially  distributed  industrial  explosion,  so  that  one  would  expect 
the  radiated  seismic  wave  to  be  depleted  in  high  frequency  spectral  content  compared 
to  that  from  a  comparable  spatially  concentrated  explosion.  Therefore,  the  directly 
radiated  compressional  (P)  wave  from  a  large  industrial  explosion  should  appear 
earthquake-like,  that  is  reduced  in  high  frequency  amplitude,  relative  to  comparable 
nuclear  tests  of  either  the  tamped  or  decoupled  types. 

Joint  Seismic,  Acoustic  and  Electromagnetic  Discrimination  Approach 

The  use  of  both  the  EM  and  seismic  methods  together  to  identify  seismic  sources 
is  summarized  in  Table  1.  Here,  in  the  column  on  the  left,  we  compare  the  expected 
signal  levels  (relative  to  noise)  for  seismic,  acoustic  and  EM  field  detections  from 
different  source  types  listed  in  the  row  along  the  top. 


Multiple  Field  Discrimination  of  Small  Seismic  Events 
A.  List  of  Expected  Qualitative  Differences _ 


Regional 

High  Frequency  (f  >  10) 
Seismic  P  wave  Signal 
(Relative  amplitude)  _ 


Small  (m^  <  4)  Event  Types,  all  with  the  same  “low"  frrqetmcy  seismic  P  ware  Amplitude 


Tamped 
Nuclear  Test 


Moderate  /  Large 


Decoupled 
Nuclear  Test 


Large 


Industrial 

Explosion 


Small  /  Moderate 


Shallow 

Earthquake 


Small 


Near  -  Regional 
Acoustic  Signal 

(Relative  amplitude) 


Small  / 
Undetectible 


Undetec  tible 


Large 


Ionospheric  E-M 
Sounding 

(Relative  amplitude) 


Small 


Undctectible 


Large 


Secondary  Ionospheric 

E-M  Emissions 

(Relative  amplitude) 


Small 


Undeicclible 


Large 


Undetectibie 


Small  / 
Undetectibie 


Small  / 
Undeieciibk 


f 

mplitude 

Scale 


Noise 

t.erel 


M 

S 

o 


signature: 
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signature: 
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signature: 
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These  event  types  are  comparable  in  that  they  are  chosen  to  have  the  same  low  fre¬ 
quency  (f  <  1Hz)  P  wave  signal  level.  (That  is  they  are  chosen  to  have  the  same  stan¬ 
dard  mb  value.)  For  such  a  set  of  "normalized"  event  types  we  list  the  signal  levels,  in 
qualitative  terms,  to  be  expected  for  each  of  the  signal  types,  in  the  table  entries. 
These  provide  an  expected  "signature"  for  each  of  the  event  types.  The  method  of 
discrimination  of  the  events  is  then  indicated  by  the  arrows  showing  what  comparisons 
could  be  made  between  "signatures"  to  identify  the  events. 

In  the  following  sections  we  describe  current  results  involving  both  observational 
studies  and  theoretical  modeling  studies  of  the  physical  phenomena  underlying  these 
discrimination  methods. 

II.  High  Frequency  Seismic  Signal  Detection  and  Spectral  Discrimination:  Obser¬ 
vations 

Spectral  discrimination  has  been  tested  in  the  past  (eg.  Savino  et.  al.,  1980)  and 
found  to  be  effective  for  events  of  magnitude  larger  than  about  3.5  to  4.0.  These  tests 
were  conducted  using  relatively  low  frequencies  (ie.  using  P  wave  magnitudes  meas¬ 
ured  at  about  .3  Hz  versus  those  at  about  3.5  Hz.)  The  method  was  termed  a  "Variable 
Frequency  Magnitude"  method  or  "VFM"  method.  Figure  1  shows  an  example  of  sin¬ 
gle  station  discrimination  in  the  far  regional  distance  range.  This  method  has  been 
proposed  by  Evemden  et.  al.  (1986)  for  the  discrimination  of  small  events,  in  the 
range  mb  ~  2  to  mb  ~  4,  at  regional  distances.  This  requires  measurement  of  high  fre¬ 
quency  data  in  the  band  from  about  1  Hz  out  to  25-30  Hz.  Based  on  preliminary 
observations  of  high  frequency  P  wave  data  in  the  near  regional  distance  range  it  was 
inferred  that  such  a  high  frequency  band  could  be  obtained  at  stations  in  stable 
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EXPLOSIONS 


-8- 


Figure  do.).  Distribution  of  events  from  the  Eurasian  data  set  in  the  rcgionaJ  distance  range  from  the  station  KAAO  at 
Kabul,  Afghanistan  (solid  triangle  on  map.)  The  solid  circles  denote  events  identified  as  nuclear  explosions  based  on  v'FM 
results  from  KAAO  and  other  stations  in  the  detection  network.  (Savino  el.  a!..  1980.) 


mb  (0.55  Hz) 


Figure  (tb).  Regional  VFM  results  obtained  at  the  station  KAAO  for  the  event  set  indicated  in  Figure  (23)  The  even 
magnitudes  range  from  somewhat  more  than  mb  =  3.5  to  somewhat  less  than  mb  =  6.0.  The  separation  of  nuclear  cxploMons 
and  earthquakes  is  of  the  order  of  .75  magnitude  units.  Nine  of  ten  explosions  were  tests  at  the  Kazakh  test  site  while  one 


was  a.  the  north  end  of  the  Casptan  sea  (No.  22).  All  the  explosions  and  earthquakes  had  similar  propagation  paths  to  the 

f  “"8  ” 7"“  •"*'  ™  «  no.  like.,  .o  be  doc  10  pad,  d.ffc, cnees  bc^cTw Tuas“, c 

ndmed  bv  die  iw„  ""S  “C'  'hC'c,‘’"-  l,kel>'  10  me-“"'cs  of  basic  special  differences  in  .he  compicss.onal  wa.es 

ramatul  by  the  two  source  types  (Savino  el  a 1 .  1980.) 
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continental  areas  out  to  about  1000  km.  and  in  tectonic  areas  to  around  300-400  km. 


Data  obtained  from  the  (former)  Soviet  Union  near  the  Kazakh  test  site  showed 
encouraging  results.  In  particular,  high  frequency  signals  from  chemical  explosions 
(m5  ~  2)  were  well  recorded  at  distances  as  large  as  640  km.  (The  maximum  distance 
of  recording.)  Figure  2  shows  an  example  recorded  at  a  station  (KSU)  with  noise 
characteristics  that  were  not  particularly  good.  (Better  results  could  be  expected  at  a 
relocated  site  in  the  same  area.)  In  any  case,  the  spectral  range  for  the  P  wave  signal 
with  S/N  >  1  was  from  about  .7  to  28  Hz. 

Data  from  other  explosions  (nuclear  and  chemical)  were  also  observed  and  some 
characteristics  were  summarized  by  Givens  et.  al.  (1989),  as  shown  in  Figure  3.  The 
data  came  from  stations  involving  both  tectonic  and  stable  continental  paths,  so  that 
there  is  scatter  in  the  distance  trends  of  the  high  frequency  "cut-off'  for  this  reason. 
In  addition,  the  yields  of  the  various  explosions  were  very  different.  Nevertheless 
there  were  indications  that  one  can  expect  to  observe  30  Hz  P  wave  signal  at  distances 
of  from  700  to  1000  km.  from  explosion  sources  as  small  as  10  tons,  and  to  20  Hz  at 
distances  of  the  order  of  about  1200  km.,  for  paths  in  stable  continental  areas. 

Our  recent  results,  relating  to  the  efficiency  and  nature  of  very  high  frequency  P 
wave  propagation  in  tectonic  and  stable  continental  areas,  suggest  that  previous  infer¬ 
ences  by  Evemden  et.  al.  (1986)  were  too  crude  and  did  not  properly  characterize  the 
high  frequency  P  wave  signal  propagation.  In  these  estimates  the  approximation  for 
the  ratio  of  P  wave  spectra  at  two  distances  A  and  Aq  is  given  by  (for  Pn  or  P): 

,4(«,A)//t„(a,,A„)  =  (A0/A)"eipi-;?4;A  -  A„)] 
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Seismograms  of  components  of  ground  motion  recorded  at  a  distance  of  637 
km.  from  a  10  ton  chemical  explosion  in  the  area  of  the  nuclear  test  site  near 
Semipalatensk  in  Kazakhstan. 
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Fig.  (2b):  Recording  of  vertical  component  of  ground  motion  from  a  10  ton  chemical 
explosion  at  a  distance  of  650  km,  in  the  area  of  the  Kasakh  test  site,  USSR. 
(From  NRDC-Soviet  Academy  Nuclear  Test  Monitoring  Program,  Sept. 
1987.) 
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Raw  (uncorrected  for  system  response)  spectra  of  4  seconds  of  signal  and 
noise  from  the  ten  ton  chemical  explosion  shown  in  Fig.(2b),  at  a  distance  of 
850  km.,  recorded  in  the  area  of  the  Kazakh  Test  Site,  USSR.  The  lines  on 
the  frequency  scale  at  the  bottom  indicate  (equally  spaced)  frequency  sam¬ 
ples  in  the  spectral  range  where  the  signal  to  noise  ratio  is  greater  than 
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Figure  (3).  Maximum  measurable  frequency  (S/S  >  2)  in  P  wave  signals  versus  distance 
for  various  sized  explosions  (From  Givens  et.  al.,  1989). 


with  A  denoting  epicentral  distance  and  Qp  and  Vp  being  the  path  averaged  dissipation 
function  and  average  compressional  velocity  respectively.  The  geometric  exponent 
factor,  m,  is  two  for  Pn  in  a  very  simple  layered  geometry,  but  in  general  can  be 
expected  to  be  in  the  range  from  about  1/2  to  3  in  complex  structures.  An  approxi¬ 
mate  form  adopted  for  Qp(f)  in  the  high  frequency  range  is: 

Qp  (0  =  Qo 

In  the  study  by  Evemden  et.  al.  1986  a  geometric  fall-off  varying  inversely  with 
the  square  of  distance  was  used  and  a  very  high  P  wave  Q,  in  the  range  near  9000, 
was  inferred  from  eastern  U.S.  data.  This  was  used  as  a  standard  for  stable  continental 
P  wave  propagation,  that  is  "Pn",  in  the  regional  distance  range.  However,  results 
from  USSR  studies  of  JVE  data  show  that  the  geometric  effect  is  closer  to  an  exponent 
dependence  of  -1  on  distance,  or  even  -.75  on  some  paths,  and  that  the  Q  for  ”Pn"  is  in 
the  range  from  to  3000  to  4000  in  stable  areas.  These  results  are  illustrated  in  Figure  4 
where  the  best  fits  to  the  amplitude  ratio  data  from  the  JVE  event  is  given  by 
geometric  spreading  exponents  (m)  near  1.  The  geometric  exponent  for  tectonic  areas 
is,  on  the  other  hand,  closer  to  2  while  the  Q  value  turns  out  to  be  in  the  range  from 
1000  to  2000,  and  typically  higher  than  the  1000  value  used  by  Evemden  et.al.  This 
is  illustrated  by  the  results  in  Figure  5. 

The  net  result  is  that  we  may  be  able  to  expect  detection  capability  of  very  high 
frequency  P  waves  out  to  the  distances  used  by  Evemden  et.  al.,  but  the  propagation 
characteristics  are  actually  somewhat  different  than  those  inferred  by  them.  Neverthe¬ 
less,  spectral  discrimination  at  regional  distances  should  be  possible  to  mb  ~  2  for 
different  kinds  of  explosions  and  earthquakes. 
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Figure  (4)  Fits  of  observed  spectral  ratios  of  PB  and  P  from  the  J VE  at  ARU  (Arti,  U.S.S.R)  and  OBN  (Obninsk,  U.S.S.R.)  normalizet 
by  the  P.  spectra  at  MSD  (Bayanaul,  U.S.S  Jt.)  near  the  Soviet  test  site.  A  better  fit  to  the  PB  ratio  at  ARU  is  obtained  fc) 
reducing  the  spreading  exponent  m  to  1.5.  At  OBN  the  data  shows  scatter,  but  the  spreading  exponent  is  clearly  less  than  1 
A  best  fit  is  obtained  for  a  mildly  frequency  dependent  Qp.  Both  paths  are  in  a  stable  continental  area  and  the  mean  Q, 
values  are  quite  high  at  high  frequencies. 


A  /  Ao  A  /  Ao 


Pn/P  Wave  Attenuation  at  Fixed  Frequency  versus  Distance 
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Pn/P  Wave  Attenuation  at  Fixed  Frequency  versus  Distance 
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Figure  5  Theoretical  fits  to  spectra  ratio  data  for  Pn  between  the  station  CHS 
Garm  and  a  temporary  station  near  the  Soviet  Kazakh  test  site. 
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We  conclude  that  on  the  basis  of  these  observations  the  predictions  of  regional 
discrimination  by  VFM  methods  should  be  effective  to  mb  «  2  at  distances  of  the  order 
of  1000  km  in  stable  continental  areas.  Station  distances  would  have  to  be  lower  in 
tectomc  areas,  around  500km,  to  achieve  such  low  magnitude  identification  capability, 
but  since  such  areas  are  limited  globally  this  does  not  imply  that  extremely  large 
numbers  of  stations  would  be  required  for  coverage.  Indeed,  the  station  distributions 
for  monitoring  a  land  mass  as  large  as  the  former  Soviet  Union  would  be  about  the 
same  as  that  estimated  by  Evemden  et.  al.  (1986). 

III.  Modeling  High  Frequency  Seismic  Signal  Propagation 

Evaluation  of  various  computational  methods  for  seismic  wave  field  propagation 
has  been  the  first  objective  of  this  part  of  our  investigations.  Useful  techniques 
include  standard  finite  difference,  finite  element,  and  pseudospectral  methods.  While 
standard  finite  difference  and  finite  element  techniques  allow  straightforward  and  rea¬ 
sonably  accurate  treatment  of  free  surface  and  absorbing  boundary  conditions  on  com¬ 
putational  grid  edges,  they  are  prohibitively  expensive  for  application  to  wave  propaga¬ 
tion  to  regional  distances  and  up  to  the  frequencies  of  interest.  A  typical  second  order 
finite  difference  method  requires  a  computational  grid  fine  enough  to  sample  the  smal¬ 
lest  wavelength  disturbance  with  8  or  10  grid  points  in  order  to  avoid  numerical 
dispersion  and  consequently  results  in  a  very  large  computer  memory  requirement  and, 
as  well,  long  computational  times. 

The  primary  advantages  of  a  pseudospectral  method  are  reduced  core  memory 
requirements,  faster  run  times,  and  improved  accuracy  for  solutions  containing  large 
gradients.  Because  the  method  may  require  only  2  or  3  computational  grid  points  per 
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smallest  wavelength  to  avoid  numerical  dispersion,  it  can  be  up  to  25  times  as  efficient 
as  a  second  order  finite  difference  solution  of  equal  accuracy  for  2-dimensional  prob¬ 
lems.  The  actual  sampling  requirements  depends  on  the  variations  of  the  material 
structure. 


The  Pseudospectral  Method 

For  our  purposes,  a  pseudospectral  method  refers  to  solution  of  a  partial 
differential  equation  by  employing  spectral  estimation  of  the  spatial  dependencies  and 
finite  differencing  in  time.  Spatial  differentiation  is  performed  in  the  spectral  domain 
and  then  transformed  back  to  the  spatial  domain.  Typically  the  field  variables’  spatial 
dependence  is  expanded  in  a  Fourier  basis,  so  that  use  can  be  made  of  efficient  Fast 
Fourier  Transforms  in  obtaining  a  spectral  representation. 


Consider  a  first-order  spatial  derivative  of  a  field  variable  p  in  one  dimension,  x: 

— p  (x,t)  is  evaluated  by  applying  a  spatial  Fourier  transform  to  p(x,t)  at  discrete 

locations  corresponding  to  computational  grid  points  in  the  x-direction.  After  multi¬ 
plying  the  transformed  quantity  p  (k^t)  by  wave-numbers  ik„,  n  =  1,...,  N,  where  N  is 


the  number  of  grid  points  in  the  x-direction,  the  derivative  is  obtained  by  performing 


the  inverse  tranform.  For  mixed  partial  derivatives,  e.g.  — p(x,y,t),  the  process  is 
repeated  along  the  other  dimension. 


By  expanding  the  dependent  variables’  spatial  dependence  in  a  trigonometric 
basis,  optimal  grid-point  per  wavelength  sampling  is  attained.  The  Nyquist  criterion 
indicates  that  theoretically  only  two  gridpoints  are  required  per  minimum  wavelength 
within  the  grid.  One  may  view  the  pseudospectral  method  as  either  the  trigonometric 
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interpolation  of  the  derivative  of  a  function  or  as  a  finite  difference  algorithm  of  max¬ 
imum  order,  where  all  neighboring  grid  points  along  one  grid  dimension  are  used  in 
the  local  computation  of  the  derivative.  The  two  interpretations  are  equivalent. 

For  a  linearly  elastic,  heterogeneous  medium,  we  formulate  the  equations  for  par¬ 
ticle  motions  as 

Pvi  =  Tij,j  +  Pfi  (1) 

and 

fij  =  +  Vjj),  (2) 

where  p.  A,  and  |!  are  the  position-dependent  density  and  elastic  constants  of  the 
medium,  and  v,  and  x,j  denote  particle  velocities  and  stresses  beyond  the  equilibrium 
configuration.  Dots  indicate  time  derivatives,  ami  commas  between  indices  denote  par¬ 
tial  differentiation  with  respect  to  the  spatial  coordinate  labeled  by  the  index  following 
the  comma.  5^  is  the  Kronecker  delta  function  and  fj  corresponds  to  body  force  com¬ 
ponents.  Equations  1  and  2  are  solved  by  a  pseudospectral  method,  and  appropriate 
boundary  conditions  must  be  applied.  Treating  the  stress  tensor  explicitly  allows  for 
incorporation  of  anelastic  rheologies  and  explicit  specification  of  a  free  surface  boun¬ 
dary  condition  on  one  side  of  the  grid.  The  periodic  nature  of  the  boundaries  of  a 
Fourier  pseudospectral  formulation  prohibit  the  use  of  absorbing  boundary  conditions 
on  the  other  sides,  so  damping  must  be  used  along  grid  boundaries  to  avoid  reflections. 

Because  the  accuracy  of  the  pseudospectral  method  surpasses  that  of  a  lower- 
order  differencing  scheme,  it  is  comparatively  more  useful  for  approximating  solutions 
that  contain  large  gradients  or  discontinuities.  Nevertheless,  large  gradients  or  discon- 
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tinuities  in  the  true  solution  cause  Nyquist  spatial  wavenumber  contamination  of  the 
pseudospectral  solution.  At  the  Nyquist  frequency  the  spectrum  of  a  real-valued  space 
variable  should  have  no  imaginary  component.  However,  the  phase  shift  introduced  by 
the  first-order  derivative  operator  ik  causes  imaginary  contributions  at  the  Nyquist  fre¬ 
quency,  effectively  violating  the  FFT’s  assumptions  of  symmetry  and  real-valued  input. 
This  error  is  particularly  serious  for  the  modelling  of  surface  and  interface  waves, 
which  propagate  along  discontinuities.  In  order  to  accurately  determine  the  effects  of 
lateral  variations  in  material  structure  on  surface  wave  phases,  it  is  important  to 
resolve  or  at  least  minimize  this  error. 

The  Chebychev  Expansion  Method 

Some  authors  have  eliminated  the  Nyquist  error  by  formulating  equations  1  and  2 
with  the  field  variables  and  material  parameters  staggered  within  each  computational 
grid  cell  (Witte  and  Richards,  1990).  In  a  pseudospectral  method,  the  relative  space 
between  the  location  at  which  variables  are  evaluated  adds  a  phase  shift  to  each  vari¬ 
able  in  the  spectral  domain,  effectively  concelling  the  phase  shift  due  to  the  derivative 
operator  and  ensuring  no  imaginary  components  at  Nyquist  frequency.  We  have  found 
such  formulations  to  be  accurate  for  body  wave  calculations,  but  surface  waves 
develop  a  well-known  ringing  as  the  solution  is  integrated  in  time.  The  staggered-grid 
formulation  appears  to  make  this  problem  worse. 

Other  authors  have  successfully  eliminated  surface  wave  ringing  by  choosing  a 
computational  grid  with  the  spacing  reduced  near  the  free  surface.  In  particular,  one 
choice  is  a  grid  whose  spacing  with  depth  corresponds  to  the  zeroes  of  Chebychev 
polynomials  of  increasing  order.  The  use  of  this  functional  form  is  desirable  in  that 
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the  Chebychev  basis  can  be  expressed  in  terms  of  trigonometric  functions,  retaining 
the  use  of  Fast  Founer  Transforms  in  the  solution  (Kosloff,  et.  al.  1990).  The  finer 
gnd  spacing  makes  available  higher  wave-numbers  for  synthesizing  the  wavefield  in 
the  vicinity  of  the  free-surface  discontinuity,  and  the  non-periodic  nature  of  the  basis 
allows  the  use  of  absorbing  boundary  conditions  on  the  opposite  end  of  the  grid. 

In  Figure  6  we  compare  the  analytic  solution  of  Lamb’s  problem  to  that  of  the 
Chebychev  method  for  a  vertical  vector  surface  source  with  a  maximum  frequency 
content  of  1  Hz.  The  material  paramenters  correspond  to  a  compressional  velocity  of 
5.4  km/s,  a  shear  velocity  of  3.1km/s,  and  a  density  of  2.5g/cm3.  The  grid  is  scaled  so 
that  the  grid  points  in  the  vicinity  of  the  free  surface  sample  between  10  and  30  times 
per  wavelength.  The  numerical  solution  compares  well  to  the  exact  analytic  result. 

The  primary  disadvantage  to  the  Chebychev  basis  method  is  that  a  stable  solution 
requires  a  fourth-order  Runge-Kutta  method  for  integration  in  time.  This  makes  the 
algorithm  slow.  We  are  presently  looking  into  other  time-integration  schemes  for  the 
Chebychev  method,  in  addition  to  alternative  solution  techniques  on  a  regular  grid. 

The  Complex  Pseudospectral  Method 

As  an  alternative  technique  for  avoiding  the  Nyquist  error,  we  have  chosen  to 
solve  the  equations  of  motion  on  a  non-staggered  grid  but  now  define  each  field  van- 
able  to  be  complex.  Solving  for  both  real  and  imaginary  components  at  each  time  step 
preserves  all  phase  information  and  circumvents  the  I  lyquist  noise  problem  in  the  same 
manner  as  does  the  staggered-gnd  technique. 

Figure  7  shows  the  Rayleigh  wave  produced  from  Lamb’s  problem  using  a 
complex-field  pseudospectral  method  and  the  same  matenal  parameters  and  source  as 
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Figure  6-a:  Lamb’s  Problem:  Comparison  of  analytic  and  Chebychev  solutions. 
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Figure  7-a:  Lamb's  Problem:  Comparison  of  analytic  and  complex  pseudospectral  solutions. 


for  the  solution  in  Figure  6.  The  grid  spacing  was  chosen  to  give  only  3  nodes  per 
wavelength  of  the  Rayleigh  wave.  The  characteristics  of  the  complex  Fourier  pseudos- 
pectral  solution  are  different  than  those  of  the  Chebychev  solution,  but  it  appears  to  be 
qualitatively  of  similar  accuracy.  However,  there  is  a  somewhat  larger  error  in  the 
Rayleigh  wave  arrival  time  (which  increases  with  increasing  distance)  than  for  the 
Chebychev  approximation. 

The  complex  formulation  requires  doubling  the  storage  for  the  velocity  and  stress 
variables,  but  the  rapid  use  of  complex  FFT’s  can  be  retained  without  additional  com¬ 
putational  effort.  In  contrast  to  the  Chebychev  method,  standard  second-order  accurate 
time  integration  is  possible  with  the  complex  Fourier  method.  This  allows  larger  time 
steps  and  a  much  simpler  and  therefore  rapid  time  integration  scheme. 

Tests  of  the  complex  pseudospectral  (CPS)  method  can  also  be  made  relative  to 
ID  modal  solutions  (eg.  Harvey,  1981).  Here  we  expect  the  numerical  CPS  result  to 
be  quite  accurate  for  body  waves  but  that  it  may  produce  some  error  in  the  fundamen¬ 
tal  mode  surface  wave  due  to  the  abrupt  velocity  change  represented  by  the  free  sur¬ 
face.  The  Modal  theory  will  of  course  produce  a  very  accurate  representation  of  the 
low  order  surface  wave  modes  so  that  this  comparison,  in  a  representative  layer  model 
of  the  earth’s  crust  and  upper  mantle,  should  give  us  a  good  check  on  the  accuracy  of 
the  CPS  method  for  "realistic"  models. 

Figure  (8)  shows  a  Modal  -  CPS  comparison  for  a  layered  earth  model  represen¬ 
tative  of  the  structure  in  the  nuclear  test  site  area  in  Kazakhstan.  The  agreement 
between  the  two  computations  for  the  surface  wave  is  close,  and  the  remainder  of  the 
seismogram  is  also  in  good  agreement.  We  therefore  consider  the  CPS  method  to  be  a 
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good  and  quite  reliable  means  of  computing  seismic  synthetics  in  complicated  media. 
Nevertheless,  we  expect  to  be  able  to  improve  on  it  and  eliminate  the  (relatively  small) 
group  delay  errors  by  one  of  several  modifications  now  under  study. 

Seismic  Modeling  Results  for  Complex  Earth  Models 

At  present  ID  analytical  methods,  employing  mode  superposition  methods,  as 
described  by  Harvey  (1981,  1991),  are  also  being  employed  to  compute  synthetic 
seismic  time  series  from  a  variety  of  source  types.  These  are  being  compared  to 
numerical  modeling  results  as  a  check  and,  more  importantly,  to  observed  data.  Here 
our  objective  is  to  systematically  vary  1  D  structures  in  a  manner  designed  to  test  the 
limits  of  ID  modeling  of  regional  seismic  data.  In  this  endeavor  we  are  using  "verti¬ 
cal  randomization"  of  the  structure  (that  is,  superimposed  random  fluctuations  on  a 
mean  velocity  model  that  fits  travel  time  data)  to  try  to  fit  near  field  and  regional 
observations  from  relatively  simple  explosion  sources. 

Figure  (9),  is  an  example  of  such  comparisons.  It  is  evident  that  by  the  use  of 
superposed  vertical  randomization  one  can  obtain  reasonable  first  order  fits  to  the  data 
over  at  least  the  0  to  1Hz  frequency  band.  However,  this  is  only  the  case  in  the 
regional  distance  range  out  to  200-300  km,  at  least  near  the  Soviet  East  Kazakh  test 
site,  and  only  over  the  lower  frequency  range  up  to  a  few  Hertz.  In  particular,  it  has 
become  apparent  that  lateral  variations  become  much  more  important  when  such  sim¬ 
ple  ID  model  results  are  systematically  compared  to  more  distant  observations  and/or 
compared  to  higher  frequency  data.  It  appears,  at  least  in  the  southwestern  Soviet 
Union,  that  lateral  variations  produce  marked  effects  on  the  signal  beyond  distances 
corresponding  to  about  fifty  wavelengths  in  propagation  distance.  (This  of  course 
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Figure  (9.) 


Comparison  of  observed  (MSD 
Verification  Explosion  (JVE)  at  1. 
crate  the  synthetic  is  shown  in  the  I 


vanes  with  wave  type,  with  the  Pg,  Lg  and  Rg  signals  understandably  showing  the 
greatest  sensitivity  to  lateral  variations.) 

Currently  we  are  developing  modal  synthesis  methods  for  2  and  3D  structures 
(Archambeau,  AFGL  annual  report,  1990)  employing  lateral  propagator  methods. 
These  will  complement  the  purely  numerical  modeling  methods  described  earlier  and 
can  result  in  marked  advantages  in  computational  speed  and  flexibility  and  may,  it  is 
hoped,  be  particularly  useful  in  inversion  procedures. 

Both  the  numerical  and  analytic  (modal)  methods  can  be  used  to  evaluate  the 
effects  of  simple  random  spatial  fluctuations  in  material  properties  in  and  vertical  and 
horizontal  directions  within  a  layered  earth  model.  Such  computations  would  be  a  first 
step  in  assessing  the  importance  of  medium  fluctuations  and  larger  scale  material  pro¬ 
perty  variability  in  the  production  of  seismic  wave  fields  seen  in  the  regional  distance 
range.  Other  more  extreme  and  complex  models  can  also  be  explored,  but  a  logical 
first  step  it  to  investigate  simple  random  variability  from  layered  mean  velocity  models 
of  the  earth. 

In  this  regard,  we  note  that  Figure  (9)  compares  observed  data  and  a  synthetic 
seismogram  generated  from  a  1-D  vertically  randomized  mean  velocity  model.  The 
mean  velocity  model  was  constrained  to  fit  the  travel  time  data  obtained  from  studies 
in  the  region  while  the  randomization  used  corresponds  to  fluctuations  in  seismic  velo¬ 
cities  and  densities  about  there  mean  values,  of  about  15%  near  the  surface  to  a  few 
percent  at  depths  in  the  upper  mantle  just  below  the  crust-mantle  boundary. 

The  fit  of  the  synthetic  to  this  single  observed  seismogram,  recorded  at  a  distance 
of  253  km.  from  the  JVE  nuclear  test,  is  certainly  not  exact  but  has  a  structure  that  is 
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quite  similar  to  that  observed.  The  lower  trace  is  an  overlay  of  the  synthetic  on  the 
observed  time  series.  The  velocity  structure  used  is  shown  in  the  lower  inset. 

As  noted  earlier  however,  test  computations  at  greater  distances,  using  the  same 
or  similar  velocity  structures,  does  not  produce  good  fits  to  observed  data  from  this 
same  event.  In  particular,  the  Rg  wave  is  observed  to  attenuate  very  rapidly  with 
increased  distance  and  the  Rg  predictions  do  not  attenuate  nearly  so  rapidly  and  the 
observed  Pg  and  Lg  are  larger  than  those  predicted.  We  conclude  from  this  that  simple 
vertical  randomization  in  the  ve.ocity  structure  is  not  sufficient  to  explain  the  basic 
first  order  properties  of  the  data  observed. 

The  next  step  is  therefore  to  include  a  form  of  lateral  variability  in  the  earth 
models.  The  simplest  variation  that  can  be  expected  is  lateral  fluctuations  in  the 
medium  velocities  as  well  as  the  vertical  fluctuations.  Examples  of  the  evaluation  of 
models  of  this  type  are  shown  in  Figures  (10),  (11)  and  (12).  These  results  were  gen¬ 
erated  using  the  CPS  numerical  method. 

The  various  model  fluctuations  indicated  in  these  figures  are  extreme  enough  to 
produce  a  range  of  results  that  allow  an  assessment  of  the  role  of  material  property 
fluctuations  on  seismic  signal  data;  particularly  the  magnitude  of  the  scattering  effects 
on  Rg  ,  Lg  ,  and  Pg. 

These  results  indicate  that  the  likely  spatial  fluctuations  in  seismic  velocities  within  the 
earth  are  not  sufficient,  by  themselves,  to  explain  the  extreme  attenuation  of  Rg  with 
distance  nor  the  large  magnitude  of  the  Lg  and  Pg  wave  trains. 

It  is  therefore  likely  that  large  scale  variations  in  shallow  structure,  coupled  with 
material  property  fluctuations  and  anelastic  attenuation  are  in  combination  responsible 
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Figure  10-a:  Synthetic  seismograms  for  a  layered  earth  with  horizontal  and  vertical  randomization 

of  1km.  square  segments  of  the  earth  structure.  The  fluctuations  very  smoothly  from 
20%  at  the  surface  to  3%  at  a  depth  of  3km. 
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Figure  10-b:  Cross-sectional  representation  of  the  randomization  applied  to  the  layered  structure. 


Figure  1 1-a:  Synthetic  seismograms  for  a  layered  earth  with  horizontal  randomization  of  4km-wide 

segments  of  the  earth  structure.  The  fluctuations  vary  smoothly  from  15%  at  the  sur¬ 
face  to  3%  at  a  depth  of  3km. 


Figure  1 1  -b:  Cross-sectional  representation  of  the  randomization  applied  to  the  layered  structure. 
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Synthetic  seismograms  for  a  layered  earth  with  12km-wide  horizontally  and  vertically 
randomized  segments  placed  every  36km.  The  fluctuations  vary  smoothly  from  15% 
at  the  surface  to  3%  at  a  depth  of  3km. 


Figure  12-b:  Cross-sectional  representation  of  the  randomization  applied  to  the  layered  structure. 


for  the  effects  observed  in  the  propagation  of  the  major  regional  signal  groups.  The 
dominent  effect  for  Rg  and  Lg  is  most  likely  to  be  strong  large  scale  lateral  variations. 
Future  modeling  studies  will  involve  systematic  investigations  of  these  structure  effects 
on  the  seismic  data  and  we  fully  expect  to  be  able  to  isolate  the  factors  of  importance 
given  our  current  modeling  capabilities. 

As  another  example  of  the  application  of  the  seismic  modeling,  Figure  (13)  shows 
seismic  time  senes  computations  that  are  designed  to  evaluate  the  effects  of  rough  sur¬ 
face  topography  near  a  buried  explosion  source.  The  effects  of  shallow  fine  scale 
layering  are  also  indicated  by  these  examples.  Clearly  the  examples  show  that  com¬ 
plexities  in  the  P  and  Pn  coda  of  the  type  observed  in  the  earth  can  be  produced  by 
topographic  features  and/or  near  source  fine  scale  layering.  In  subsequent  applications 
a  generalized  version  of  the  current  computer  program  will  include  non-linear  effects, 
including  failure,  so  that  the  effects  of  spallation  (separation  of  near  surface  layers 
under  explosive  shock  loading)  can  also  be  included.  The  effects  of  strong  or 
moderate  lateral  variations  in  near  source  shallow  structures  can,  of  course,  also  be 
included  now  and  may  be  systematically  studied.  Through  this  approach  we  plan  to 
"sort  out"  the  quantitative  nature  of  each  of  these  effects  on  the  directly  radiated 
seismic  field  from  explosions  and  be  able  to  evaluate  their  total  impact  on  source  depth 
estimates,  yield  estimation  and  discrimination. 

In  addition  to  near  source  effects  modeling,  we  are  using  both  analytical  and  new 
numerical  methods  to  extend  our  computational  range  to  greater  regional  and  telese- 
lsmic  distances.  Most  of  these  methods  involve  2  and  3D  modeling  in  order  to  be  able 
to  account  for  lateral  variability  in  the  medium.  Currently  restrictions  imposed  by 
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figure  (13-a.)  Contours  of  the  vertical  particle  velocity  from  an  explosion  3(H)  motets  below, 
the  ftee  stnface  with  high  relief  surface  topography  in  the  vicinity  ot  the  souicc 
aiea.  The  dimensions  of  the  cross-section  shown  is  km  wide  h\  ' -C  km 
deep  Maximum  sm face  elevation  is  31  in  motets. 


(  1  lioui-.  :  T ’ «i  vt-ci..  i!  niiiiu'k'  velocity  hum  an  explosion  3(K)  meters  below 
the  lire  vn,  .1  v  ;ih  nil’ll  leliet  sintaee  topogiaphy  in  the  vicinity  of  the  souicc 
I  he  uimn'Miuv.  o»  the  eioss-scction  shown  is  «)  km.  wide  by  3.45  km 
dei  p  M.iviiiimn  -nM.K.e  elevation  is  300  meters 
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computer  memory  and  speed  limitations  confine  the  numerical  seismic  modeling  appli¬ 
cations,  in  2D,  to  distances  of  less  than  1000  km  and  to  frequencies  less  than  1-2  Hz. 
Modeling  in  2D  to  frequencies  of  20-30  Hz  is  usually  limited  to  distances  of  about  100 
km  from  the  source.  By  developing  moving  grid  methods  and  FFT  evaluations  of  par¬ 
tial  derivatives,  which  have  already  been  tested,  we  expect  to  be  able  to  model  high 
frequency  wave  fields  throughout  the  regional  distance  range  to  a  few  thousand  kilom¬ 
eters  and  the  primary  teleseismic  signals  (body  waves)  to  all  distances. 

IV.  Modeling  Atmospheric  Wave  Fields  and  Ionospheric  Electron  Density  Varia¬ 
tions  Due  to  Near  Surface  Seismic  Sources 

Because  of  the  exponential  decrease  of  atmospheric  density  with  height,  buoyant 
pulsed  gravity  waves  generated  by  surface  or  subsurface  seismic  sources  can  be  of 
appreciable  amplitude  throughout  the  atmosphere.  Furthermore,  above  100km  in 
height,  these  flow  transients  affect  the  ionospheric  E-M  fields  through  changes  in  the 
distribution  of  the  charged  particles.  The  basic  equations  governing  motions  of  the 
neutral  atmosphere  are  the  conservation  laws  of  mass,  momentum  and  energy  together 
with  the  ideal  gas  equation  of  state.  The  specific  nonlinear  continuum  equations  incor¬ 
porate  nonlinear  advective  terms  as  well  as  the  gravitational  field,  gas  compressibility, 
viscosity  effects  and  thermal  conductivity.  For  electron  motions  in  the  ionosphere  a 
first-order  continuity  equation  can  be  used  as  a  first  approximation  which  assumes  that 
electrons  move  with  the  neutral  atmosphere. 

For  modeling  purposes,  the  set  of  partial  differential  equations  for  the  atmosphere 
are  converted  to  a  corresponding  set  of  finite  difference  equations  in  order  to  effect 
numencal  integration  in  time  and  space.  The  non-linear  terms  are  treated  non-locally 
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on  the  lattice  for  stability,  effectively  controlling,  internally,  the  instabilities.  In  addi¬ 
tion,  random  velocities  and  pressures  are  attributed  to  the  inherent  fine  scale  turbulence 
in  the  atmosphere  and  are  incorporated  in  the  modeling,  as  are  mean  drift  particle 
velocities.  In  particular,  in  order  to  account  for  the  inherent  turbulence  in  the  atmo¬ 
sphere,  the  flow  variables  at  a  point  are  decomposed  into  a  mean  flow,  governing 
winds  and  a  perturbed  flow  that  incorporates  the  turbulence.  A  new  approach, 
designed  to  include  turbulence,  has  been  developed  using  random  perturbations, 
obtained  from  a  random  number  generator,  which  are  inserted  directly  into  the  finite 
difference  equations.  Turbulence  is  also  produced  by  a  random  distribution  of  tem¬ 
perature  at  the  surface  which  produce  thermal  structures  with  upward  and  downward 
flows.  Horizontal  winds,  impacting  on  a  variable  and  random  topography,  also  pro¬ 
duce  upward  and  downward  motions  which  have  a  random  stochastic  character. 

The  set  of  non-linear  partial  differential  equations  are  converted  to  a  correspond¬ 
ing  set  of  finite  difference  equations  for  numerical  integration  in  time  and  space  and 
upwind  differencing  is  used  for  first  order  spatial  gradients  with  the  advection  velocity 
terms  acting  at  the  upwind  point.  However,  if  the  velocity  operates  on  its  own  velo¬ 
city  gradient,  such  non-linear  terms  are  treated  non-locally  on  the  lattice  for  stability, 
effectively  controlling  internally  any  unstable  growth. 

There  are  at  least  three  types  of  boundary  important  to  the  modeling  of  fluid 
flows.  Specifically,  the  air-ground  surface  is  topographically  complex  with  a  turbulent 
boundary  layer  of  the  order  of  a  few  meters  at  the  interface.  At  this  boundary,  vertical 
velocities  are  random  both  in  time  and  spatial  locations.  Because  of  the  presence  of 
the  lower  boundary  layer  above  a  complex  topography,  horizontal  velocities  are  not 
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taken  as  zero  but  incorporate  winds  and  turbulence  effects.  The  top  atmospheric  boun¬ 
dary  is  open,  with  decreasing  density.  The  topmost  boundary  should  mimic  the  condi¬ 
tions  for  an  open  atmosphere  with  specific  considerations  for  buoyancy  and  field  gra¬ 
dients.  We  have  examined  various  options  including  fixing  velocities  and  densities 
and  their  gradients.  However,  we  have  adopted  the  general  open  flow  boundary, 
which  we  also  use  for  the  artificial  side  boundaries.  The  side  boundaries  are  artificial, 
due  to  grid  restrictions,  and  must  mimic  open  boundaries  that  allow  free  flow  in  either 
direction.  We  have  adopted  the  more  usual  approach,  wherein  the  dependent  variables 
are  constrained  to  stay  constant  at  these  open  boundaries. 

Explosive  sources  at  and  below  the  ground  are  simulated  by  stress  distributions 
and  velocities  along  the  earth’s  surface  and  their  resultant  effects  on  the  atmosphere 
are  integrated  upward  and  outward.  Various  velocity  sources  are  used  at  the  lower 
boundary  with  diffenng  time,  amplitude  and  radial  dependences.  The  standard  input  is 
a  source,  comprised  of  the  first  differential  of  a  gaussian  in  time,  that  approximates  the 
initial  pulse  from  an  underground  explosion.  Cartesian  coordinates  are  used  to  model 
the  3-dimensional  system,  with  the  source  at  the  center  of  the  bottom  plane. 

Conservation  Laws 

The  continuity  equations  are  based  on  the  values  of  the  fields  at  particular  points 
in  space  and  time.  Conservation  of  mass  is  expressed  as, 

(l') 

where  p  is  density  and  Uj  is  velocity  in  the  Xj  direction.  Conservation  of  momentum  is 
similarly  expressed  as: 
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where  X,  are  external  forces,  Py  is  the  generalized  stress  such  that: 


P,j  =  -  pAj  +  2peij  -  2/311-6^ 


(3.) 


with  e^  =  1/2 


3uj 

3xj 


+ 


3uj 

9xj 


the  strain  rate  and  |i  the  viscosity.  Conservation  of  energy 


is: 


(4.) 


where  T  is  temperature,  cv  is  specific  heat  at  constant  volume,  K  is  thermal  conduc¬ 
tivity,  and  4>  is  the  viscous  dissipation.  Here: 

<!>  =  2\ic?}  -  2/3  p(eij)2 

The  equation  of  state  for  the  atmospheric  gas  is  taken  to  be  ideal,  i.e. 

P  =  ^PT  (5) 

where  kB  is  Boltzmann’s  constant  and  m  is  the  mean  molecular  weight. 

Normalized  Equations 

In  the  fundamental  equations,  (1)  through  (5),  the  dependent  and  independent 
variables  can  be  normalized  with  respect  to  typical  values.  For  an  ambient  atmo¬ 
sphere,  with  exponential  decay  of  density  with  height,  distances  are  normalized 
through  the  scale  height,  H,  which,  at  the  surface,  is  approximately  8400  metres. 
Velocities  are  normalized  with  respect  to  cs,  the  sound  velocity  of  air  at  the  earth’s 
surface.  Similarly  density,  pressure  and  temperature  are  normalized  to  surface  values, 
and  the  independent  variable,  time  t,  is  normalized  by  (H/c,).  Thus  for  the  continuity 


of  mass,  we  get,  as  before. 


f  +  A(puJ,  =  0  (6.) 

where  the  new  variables  arc  now  normalized  and  given  the  same  symbol  as  the  origi¬ 
nal  variables.  Incorporating  gravity  as  the  external  force,  the  momentum  conservation 
equation,  (2),  becomes 

+  Puj"  g~  =  -G*'8(z)P'Cx  -  ^2'g^r  +  ^ ^ 


where  Gs  =  (gsH/c2)  is  a  measure  of  the  ratio  of  potential  energy  to  thermal  energy. 
A2  =  Ps/(pscs2)  is  a  measure  of  the  ratio  of  stress  energy  to  thermal  energy. 
A4  =  (is/(ps  cs)  is  the  ratio  of  viscous  to  thermal  energy, 
y  is  the  normalized  viscosity  drag  force. 


where  |I  is  normalized  to  |is,  the  viscosity  at  the  surface.  In  the  atmosphere,  (J.  is  usu¬ 
ally  taken  to  be  constant  for  the  molecular  viscosity.  In  the  case  of  conservation  of 
energy,  the  normalized  equation  is 


where 


A5  =  IV 
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m*Cv 
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A*  =  A*  K 
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where  K  is,  as  usual,  taken  constant  for  the  atmosphere.  The  equation  of  state,  on  nor¬ 
malization,  is 

p  =  p-T-/m(z)  (9.) 

kg 

where  the  ideal  equation  of  state  at  the  surface  is  ps  =  — psTg  with  ms  the  surface 

ms 

value  of  mean  molecular  weight  (29.0)  and  m(z)  the  height-dependent  normalized 
value. 

Ionospheric  Motions 

The  basic  conservation  law  of  charged  particles,  assuming  no  creation  or  annihila¬ 
tion,  is: 

3N« _ ftNa-Ujo)  (10.) 

~5C  3xj 


where  Na  is  the  number  of  particles  of  type  a  and  ujct  is  their  velocity  in  the  j’th 
direction.  The  initial  concentration  of  the  charged  particles  is  taken  to  be  time- 
independent,  with  only  a  vertical  functional  dependence,  N(z),  where  the  subscript  a 
has  been  dropped  for  the  type  of  particle.  Assuming  only  small  changes  in  this  con¬ 
centration,  the  dependence  can  be  found  from  integrating  eqn.  (10)  over  the  range  ^  to 
the  present.  To  zeroth  order,  this  concentration  change  becomes: 

6N(z,t)  =  -ju^Odt'  -  N(z)- j-^p-dt  (11.) 

The  first  term  in  (1 1)  is  the  concentration  change  due  to  the  displacement  of  the  ionos¬ 
pheric  layer,  while  the  second  term  arises  as  a  result  of  compression  or  rarefaction  and 
is  the  predominant  term  when  dealing  with  processes  involving  characteristic  dimen¬ 
sions  smaller  than  the  width  of  the  layer.  The  velocity  of  the  charged  particle  is 
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usually  assumed  to  be  identical  with  that  of  the  neutral  gas  to  zeroth  order  and  this  is 
the  velocity  that  is  used  in  the  finite  difference  calculations  for  electron  density 
changes. 

The  initial  concentration  of  electrons  is  taken  to  be  that  of  a  Chapman  distribution 
which  has  a  maximum  density  at  345  km  height  and  decreases  rapidly  below  about  90 
km  with  the  functional  dependence  on  height  defined  by: 

N(z)  =  Ncexp(  -  j  ( 1  -^-e-*))  ( 1 2.) 

Where  =  (z-hJ/H,  h<.  =  345  km,  H  =  65  km  and  Nc  is  the  normalizing  value. 

Finite  Difference  Scheme 

The  set  of  non-linear  partial  differential  equations  are  converted  to  a  correspond¬ 
ing  set  of  finite  difference  equations  for  explicit  computer  integration  in  time  and 
space.  Upwind  differencing  is  used  for  first  order  spatial  gradients,  with  the  advection 
velocity  terms  acting  at  the  upwind  point.  However,  if  the  velocity  operates  on  its 
own  velocity  gradient,  such  non-linear  terms  are  treated  non-locally  on  the  lattice  for 
stability,  effectively  controlling  internally  any  unstable  velocity  growth. 

The  updated  variable  is  projected  not  from  just  the  old  dependent  variable,  a  pro¬ 
cess  that  is  inherently  unstable,  but  from  a  distributed  smoothed  average  of  the  vari¬ 
able  at  locations  surrounding  the  specific  spatial  location.  Such  a  smoothing  method 
brings  stability  to  the  differencing  scheme.  However,  the  attendant  numerical  diffusion 
is  minimized  by  not  smoothing  the  density  variable,  which  has  only  a  small  effect  on 
stability.  This  approach  also  helps  in  stabilizing  the  integration  at  grid  comers  and 
boundaries.  The  second  order  derivatives  in  the  viscosity  and  thermal  conductivity 
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terms  are  modeled  by  finite  differences  taken  at  the  surrounding  spatial  locations. 

In  the  explicit  integration  scheme,  the  updated  flow  velocities,  temperature  and 
density  are  obtained  via  their  continuity  equations  while  pressure  is  obtained  from 
insertion  of  the  updated  density  and  temperature  into  the  ideal  gas  equation. 

Boundary  Conditions 

There  are  at  least  three  types  of  boundaries  important  to  the  modeling  of  fluid 
flows.  The  air-ground  surface  is  topographically  complex  with  a  turbulent  boundary 
layer  attached.  The  top  atmospheric  boundary  is  open  to  space  with  decreasing  den¬ 
sity.  The  side  boundaries  are  artificial,  due  to  grid  restrictions,  and  mustjmimic  open 
boundaries  that  allow  free  flow  in  either  direction. 

At  the  bottom  boundary  vertical  velocity  functions  are  input  as  sources  of 
momenta  at  spatial  locations.  Otherwise,  as  usual,  vertical  velocities  are  taken  to  be 
zero  at  the  bottom.  Because  of  the  presence  of  the  lower  boundary  layer  above  a  com¬ 
plex  topography,  horizontal  velocities  are  not  taken  as  zero  but,  instead,  constant  velo¬ 
city  and  density  gradients  are  assumed  in  the  vertical  direciton.  For  subsurface  sources 
only  momentum  inputs  are  considered  at  this  bottom  boundary.  For  sources  at  or 
above  the  surface,  both  momentum  and  pressure  conditions  must  be  applied,  just  as  in 
atmospheric  sources. 

The  top-most  boundary  should  mimic  the  conditions  for  an  open  atmosphere  with 
specific  considerations  for  buoyancy  and  field  gradients.  We  have  examined  various 
options  including  fixing  velocities  and  densities  and  their  gradients.  However,  we  have 
adopted  the  general  open  flow  boundary,  much  as  we  use  for  the  artificial  side  boun¬ 
daries.  As  usual,  in  order  to  preserve  conservation  relations,  all  normal  gradients  are 
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set  to  zero  at  these  open  boundaries.  However,  this  would  not  permit  heat  flow 
through  the  boundary.  Therefore  the  second-order  normal  derivative  of  temperature  is 
made  constant. 

Turbulence  Effects 

In  order  to  account  for  the  inherent  turbulence  in  the  atmosphere  below  the  ther¬ 
mopause,  the  flow  variables  at  a  point  can  be  decomposed  into  a  mean  flow  and  a  per¬ 
turbed  turbulent  flow.  In  the  momentum  equation,  additional  components  are  thus 
obtained  for  the  generalized  stress.  These  are  mainly  interaction  terms  between  the 
mean  and  perturbed  densities  and  velocities,  termed  the  Reynold’s  stresses,  which 
represent  the  interaction  of  the  mean  flow  with  the  background  turbulence.  These 
extra  stresses  have  been  approximated  by  various  phenomenological  approaches. 
Boussinesq  introduced  the  concept  of  eddy  viscosities  in  order  to  use  the  Newtonian 
equations  with  the  usual  but  much  larger  viscosity  term.  Our  models  evaluate  the 
efficacy  of  this  method  using  an  eddy  thermal  conductivity.  We  also  attempt  to  evalu¬ 
ate  different  forms  of  these  Reynold’s  stresses  through  algorithmic  modeling  of  various 
drag  forces  that  mimic  the  effect  of  these  interaction  terms.  It  is  found  that  even  small 
drag  forces  of  a  particular  type  can  alter  die  flows  and  their  temporal  dependence.  An 
alternative  approach  to  turbulence  is  developed  in  the  use  of  random  perturbations, 
obtained  from  a  random  number  generator,  and  input  directly  into  the  finite  difference 
equations. 

Modeling  Results 

The  result  of  the  atmospheric  modeling  for  effects  of  a  surface  explosion  can  be  sum¬ 
marized  as  follows: 
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(I.)  A  time  dependent  transient  pulse  propagates  upward  from  the  seismic  wave 
perturbation  of  the  ground  surface  and  has  increasing  amplitude  relative  to  the 
ambient  pressure.  This  pioduces  asymmetric  flows  which  control  the  flow 
development  and  the  upward  propogation  of  the  transient.  The  initial  positive 
density  pulse  is  propagated  upward  more  slowly  than  the  following  negative  den¬ 
sity  pulse  which  has  increased  buoyancy.  This  initiates  a  sequence  of  circulation 
patterns  that  develop  through  what  appears  to  be  asymmetric  triangular  modes 
across  the  horizontal  cross-section.  The  circulation  patterns  for  the  phenomena 
are  characterized  by  upward  central  motions  of  the  lighter  matter,  which,  at  the 
neutral  buoyancy  level,  push  outward  to  the  side.  The  centroid  of  the  transient 
pulse  initially  moves  upward  rapidly,  but  slows  down  to  the  group  velocity  speed 
of  sound  in  the  atmosphere.  The  advected  air  mass  tries  to  remain  in  its  horizon¬ 
tal  stratification  in  order  to  minimize  changes  in  its  gravitational  potential.  How¬ 
ever,  it  appears  that  energy  and  momenta  are  transported  through  traveling  waves 
in  the  circulation  pattern.  Similar  effects  have  been  observed  in  the  real  atmo¬ 
sphere  when  thermals  propagate  upward  from  the  Earth’s  surface  with  similar  cir¬ 
culation  patterns. 

(2.)  After  a  model-dependent  characteristic  time  a  bifurcation  of  the  flow  occurs 
with  the  eventual  reversal  of  the  velocity  direction.  The  bifurcation  phenomena 
occurs,  in  this  model,  every  100  seconds,  so  that  it  has  a  period  of  just  over  3 
minutes.  A  drag  force  is  input  in  order  to  model  the  effect  of  the  inherent  back¬ 
ground  turbulence  of  the  atmosphere.  Such  a  drag  force,  which  removes  2%  of 
the  component  velocities  at  each  computational  grid  pount  at  each  time  step, 
removes  the  periodic  bifurcation  and  a  standing  wave  is  formed  in  the  atmosphere 
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with  constant  field  patterns.  However,  with  a  1%  removal  rate,  the  patterns  are 
periodic  with  similar  bifucations  as  in  the  zero  drag  case.  Because  existing 
atmospheric  turbulence  acts  on  the  transient  gravity  wave  as  a  perturbation,  we 
have  also  modeled  its  effect  by  imposing  a  random  component  on  each  field  at 
each  time  step  and  grid  point.  The  usual  bifurcations  are  obtained  but  with 
differing  patterns  from  the  zero  turbulence  case.  However,  the  appearance  of  the 
pressure  and  density  fields  is  more  realistic  due  to  added  diffusion  and  random 
components. 

As  the  transient  pulse  moves  upward  in  the  atmosphere,  it  magnifies  in  amplitude 
relative  to  the  exponentially  decreasing  ambient  pressure.  Thus,  the  level  at  which  a 
specific  pressure  is  located  will  oscillate  as  the  transient  pressure  pulse  moves  through. 

To  the  first  order,  the  electrons  in  the  ionosphere  are  assumed  to  move  with  the 
flow  of  the  dominant  neutrals.  Thus  the  change  in  the  electron  density  can  be  calcu¬ 
lated  from  a  conservation  law,  whose  integration  in  time  gives  the  total  electron  den¬ 
sity  variation.  The  ambient  electron  density  is  approximated  by  the  Chapman  function 
which  has  a  maximum  electron  density  at  350km  and  effectively  zero  electron  density 
below  about  90  km.  For  reasonable  synthetic  velocity  sources  at  the  ground  surface, 
we  find  that  changes  in  electron  density  from  100km  and  up  are  of  the  same  order  as 
those  observed  in  E-M  experiments  conducted  for  surface  and  subsurface  explosions. 

In  this  regard.  Figure  (14)  shows  an  example  of  the  predicted  fluctuations  in  tem¬ 
perature  and  electron  density  in  the  ionosphere  due  to  a  near  surface  underground 
explosion.  In  this  case  the  explosion  was  taken  to  be  a  tamped  underground  nuclear 
test  at  a  depth  of  300  meters  with  a  seismic  body  wave  magnitude  near  5.  (Much 
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(a.)  Time  and  spatial  variation  of  the  temperature  change  expressed  as  a  fractional  change  in  the 
temperature. 


(b.)  Time  variation  of  the  electron  density  change  expressed  as  a  fractional  change  in  the  electron 
density  at  200  km.  altitude. 


Figure  (14.)  Fluctuations  in  the  ionospheric  temperature  and  electron  density  due  to  a  transient 
pressure  pulse  at  the  earth’s  surface  from  a  near  surface  explosion. 
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smaller  industrial  explosions,  very  near  or  at  the  earth’s  surface,  would  typically  pro¬ 
duce  comparable,  or  even  larger  signals.) 

V.  Modeling  High  Frequency  Seismic  Noise:  Atmospheric  Sources 

The  nature  of  high  frequency  seismic  noise  is  indicated  by  the  observed  noise 
acceleration  power  shown  in  Figure  (15).  The  station  shown  (BAY)  is  near  the  former 
Soviet  test  site  in  Kazakh  and  is  typical  of  the  high  frequency  noise  seen  both  within 
the  former  USSR  and  elsewhere.  Three  components  of  ground  acceleration  on  the  sur¬ 
face  and  at  about  100  meters  depth  are  shown.  Both  high  and  low  wind  level  seismic 
noise  spectra  are  shown  on  each  plot.  In  all  cases  the  high  frequency  seismic  noise 
increases  with  high  wind  levels.  It  is  also  apparent  that  the  acceleration  power  is 
roughly  constant  over  the  band  from  about  1  Hz  to  30  Hz.  Above  30  Hz  the  noise 
acceleration  power  decreases  with  increasing  frequency,  particularly  in  the  bore-hole  at 
100m.  depth. 

Given  the  rather  strong  dependence  of  the  noise  level  on  wind  velocity,  it  is 
natural  to  infer  that  atmospheric  coupling  at  the  earths  surface  is  an  important  means 
of  excitiation  of  high  frequency  seismic  noise.  A  more  detailed  understanding  of  the 
atmospheric  excitation  of  seismic  noise  is  clearly  important  since  the  reduction  or 
cancelation  of  this  noise  is  dependent  on  an  understanding  of  its  origins,  mode  of  exci¬ 
tation  and  propagation  within  the  medium. 

In  order  to  investigate  the  production  of  seismic  noise  by  atmospheric  processes, 
the  atmospheric  modeling  programs  were  linked  with  the  linear  elastic  seismic  model¬ 
ing  programs.  The  lower  atmosphere,  composed  of  a  day-time  turbulence  boundary 
layer  with  a  heigh  of  2  km,  is  simulated  with  a  random  surface  topography.  Winds, 
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blowing  on  the  topography,  induce  upward  and  downward  flow  velocities.  Random 
temperature  changes  in  space  on  the  ground  surface  also  produce  flows  that  self- 
organize  into  plumes  that  coalesce  above  the  boundary  layer  into  larger  scale  thermals. 
Together  with  random  turbulence  in  the  boundary  layer,  these  flows  induce  pressure 
and  velocity  fluctuations  along  the  ground  surface.  These  effects  are  the  input  into  the 
seismic  modeling  code  which  integrates  in  time  from  the  top-most  surface  boundary. 

Preliminary  results  indicate  that  the  seismic  noise  that  is  produced  decreases  in 
amplitude  with  depth  and,  as  shown  in  Figure  (16),  produces  a  seismic  velocity  spec¬ 
trum  that  has  a  trend  that  decreases  as  1/f  with  increasing  frequency,  in  the  range  from 
about  1  to  50Hz.  Below  about  40  meters  the  seismic  noise  appears  to  interact  in  such 
a  manner  that  much  smoother  variations  in  spatial  distributions  are  obtained  than  at  the 
surface  and  with  associated  decreasing  fluctuations  in  time.  Both  topography  and 
winds  are  found  to  be  of  major  importance  in  terms  of  amplitude  and  character  of  the 
noise.  From  preliminary  results  it  can  be  expected  that  time  of  day  will  also  be  impor¬ 
tant  due  to  the  change  of  the  turbulent  boundaiy  layer  with  the  heating  of  the  Sun  and 
its  temporal  dependence. 

VI.  Summary  and  Conclusions 

We  have  tested  and  modified  our  atmospheric  modeling  capabilities  to  include  the 
most  important  non-linear  effects,  in  particular  the  effects  of  sub-grid  scale  turbulence. 
We  have  made  good  progress  on  the  turbulence  phenomena  by  introducing  a  randomly 
fluctuating  component  to  the  field  variables  (i.e.,  pressure,  velocites,  temperature  and 
density)  that  simulates  sub-grid  level  turbulent  effects.  Results  are  encouraging  and  in 
particular  give  stable  dynamical  solutions  to  test  problems  that  are  comparable 


log  amplitude 


(a.)  Seismic  noise  (velocity  versus  time)  near  the  free  surface  due  to  the  atmospheric 
turbulence. 


log  frequency  (Hz; 

(b)  Fourier  spectrum  of  the  seismic  (velocity)  noise  in  (a) .  The  trend  of  the 
mean  spectrum  is  as  f~l,  which  is  what  is  commonly  observed. 


Figure  (16).  Theoretically  predicted  seismic  noise  due  to  models  of  atmospheric 
turbulence  at  the  earth's  free  surface. 


to  observations. 


The  seismic  investigations  have  focused  on  generation  and  testing  of  2  and  3D 
finite  difference  programs  that  incorporate  surface  topography,  medium  randomness 
and  lateral  variability.  Adaptation  of  FFT  methods  coupled  with  moving  grids  have 
been  successful  when  tested  against  analytical  and  conventional  finite  difference 
methods  and  can,  in  principle,  provide  capabilities  for  predictions  of  wave  fields  in 
heterogenous  media  at  regional  distance  ranges  using  moderate  sized  computers  (e.g., 
high  level  work-stations  such  as  the  Stellar  Computer)  with  only  modest  core  size 
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requirements  (i.e.,  a  few  hundred  megabytes).  Transmitting  grid  boundaries  have  also 
been  developed  and  tested  with  success.  Analytical  (modal)  theory  methods  are  being 
developed  for  3D  laterally  varying  media,  to  be  used  along  with  the  2D  theory  already 
developed. 

Results  of  modeling  studies  and  their  comparisons  with  observed  data  have  shown 

that: 

(1)  The  atmospheric-ionospheric  modeling  predictions  of  electron  density  fluctuations 
have  amplitudes  and  wave  forms  of  the  size  and  type  inferred  from  observations. 

(2)  Coupled  atmospheric-seismic  modeling  indicates  that  atmospheric  turbulence, 
simulated  by  random  fluctuations  in  state  variables  near  the  free  surface,  produces 
high  frequency  seismic  noise  with  spectral  character  close  to  that  observed;  that  is 
with  a  velocity  amplitude  spectrum  varying  as  1/f  as  a  function  of  frequency 
above  1Hz. 

(3)  Seismic  wave  field  modeling  in  complex  structural  models,  incorporating  rough, 
near  source  topography  and  fine  scale  randomized  layering,  produces  seismic  syn- 
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thetics  having  the  complex  character  of  observed  seismograms.  Simple  vertical 
and  horizontal  randomization  can  only  be  an  adequate  representation  of  earth 
structure,  and  the  seismic  wave  fields  observed,  within  a  few  tens  of  wavelengths 
from  the  source.  Beyond  that  distance  the  effects  of  strong  shallow  lateral  varia¬ 
tions  in  both  average  and  random  characteristics  of  the  earth’s  structure  produce 
effects  in  the  observed  wave  field  that  become  of  first  order  and  therefore  impor¬ 
tant,  so  that  accounting  for  large  scale  lateral  variations  is  necessary  to  explain 
observed  seismic  wave  fields  in  the  regional  distance  range. 

(4)  Preliminary  studies  of  particular  complex  seismic  wave  types,  such  as  Lg  and  Rg, 
indicate  that  only  large  scale  lateral  variations  in  structure,  in  combination  with 
both  vertical  and  lateral  randomization  can  explain  the  wave  forms  and  attenua¬ 
tion  characteristics  observed. 

Observational  studies  that  have  been  a  part  of  our  investigations  indicate  that: 

(5)  High  frequency  seismic  P  waves  (5-30Hz)  propagate  to  considerable  distances 
with  high  efficiency  in  stable  continental  areas.  Direct  observations  indicate  that 
the  20-30Hz  band  can  be  detected  well  above  noise  levels  in  this  frequency  range 
to  distances  of  the  order  of  1000km.  In  tectonic  regions  the  efficiency  is  lower, 
with  signal  energy  in  the  20-30Hz  band  generally  detectable  to  about  400-500km. 

(6)  Quantitative  investigations  show  that  the  geometric  distance  exponent  is  near  -1 
in  stable  continental  areas  and  near  -2  in  tectonic  areas.  The  dissipation  function, 
Q,  when  constrained  to  be  frequency  independent  in  the  l-30Hz  band,  provides  a 
fit  to  P  wave  spectral  ratio  data  in  both  stable  and  tectonic  continental  areas  of 
the  former  Soviet  Union.  We  obtain  values  of  between  3000  and  4000  for  stable 
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areas  and  between  1000  to  2000  in  tectonic  areas. 


In  the  future  we  plan  to  apply  the  modeling  techniques  developed  in  a  wide  range 
of  more  systematic  studies  and  comparisons  with  observations,  as  well  as  to  extend  our 
domain  of  modeling  capability  to  greater  distance  ranges  and  higher  frequencies.  Our 
ultimate  objective  will  be  to  understand  and  predict  phenomena  under  a  variety  of  con¬ 
ditions  involving  medium  structure  and  rheological  behavior  and  to  then  devise  source 
identification  methods  based  on  these  results  while  testing  them  against  both  synthetic 
and  observed  data. 
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